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i.  nmooDcnoN 


Several  papers  in  recent  years  [1,2,3]  have  dealt  with  the  construction 
and  Migration  of  plane-wave  seismograms.  A  plane-wave  seismogram  can  be 
constructed  by  performing  a  slant-stack  over  source  or  receiver  coordinates 
with  a  constant  step-ont  in  time  per  trace.  Based  on  travel-tiaie 
considerations,  it  is  then  possible  to  conceive  of  an  iaaging  procedure 
based  on  the  concept  of  downward  continnation.  This  process  images 
structure  by  propagating  wavefields  in  reverse  back  to  the  reflectors  froa 
which  they  originated. 

It  is  not  possible,  froa  travel-tiae  considerations  only,  to  arrive  at 
the  correct  amplitude  or  weighting  factors  to  be  need  in  the  iaaging 
process.  Ideally,  the  amplitudes  on  a  migrated  or  imaged  section  shonld  be 
proportional  only  to  the  local  reflection  strength.  The  imaged  section 
shonld  be  free  froa  source-receiver  and  geometrical  spreading  factors,  and 
also  free  from  focusing  effects  due  to  the  curvature  of  the  reflecting 
surface. 

In  a  recent  paper  in  Geophysics  [3],  Trietel,  et  al,  examine  the 
question  of  amplitude  dependence  on  plane-wave  seismograms.  Their  work  was 
based  on  a  plane-wave  decomposition  technique  described  by  MBller[4].  The 
analysis  contained  in  the  Trietel  paper  is  strictly  valid  only  for 
horizontal  reflectors  located  in  a  vertically  stratified  medium.  Their 
numerical  examples,  however,  show  that  the  method  has  a  much  wider  range  of 
applicability  to  more  complex  geometries. 


This  report  outlines  the  construction  of  s  more  exset  operator  for  the 
■igration  of  plane-wave  seismograms.  The  objective  here  is  to  construct  a 
procedure  which  will  produce  an  output  appropriate  for  both  strnctnral 
inaging  and  velocity  inversion!  that  is,  an  ontpnt  which  is  a  direct 
indicator  of  reflection  strength.  The  calculations  presented  here  allow  the 
reflector  to  have  an  arbitrary  complex  two-dimensional  shape.  While  the 
reflector  is  assumed  to  be  two-dimensional,  the  propagating  wavefields  are 
assumed  to  be  three  dimensional  and  are  responses  from  point  sources.  The 
analysis  presented  here  deals  with  the  case  of  a  constant  background 
velocity  field.  More  general  cases  are  now  being  examined. 


Figure  (1)  shows  the  geoaetriesl  qosntities  involved  in  this 
discussion.  A  source  of  acoustiesl  waves  is  located  on  the  surface  at 
Xf  -  (xt,0,0),  giving  rise  to  an  incident  field  which  is  scattered  froa  the 
reflection  interface  I.  These  scattered  waves  are  then  detected  on  the 
surface  by  a  series  of  receivers  located  at  -  (xg,0,0).  The  acattered 
observations  are  denoted  by  Pg(x#,Xg,t). 

By  an  eleaentary  application  of  Huygen's  principle,  an  incident  plane 
wave  front  aay  be  synthesized  by  use  of  an  array  of  spherical  sources  along 
the  surface  which  are  separated  in  tiae  by  an  appropriate  delay.  This 
generates  a  sequence  of  point  source  responses  which  are  then  sunned  to 
synthesize  a  plane  wave  front.  In  the  tiae  donain,  this  suaaation  takes  the 
form 

Pp^(Xg,a,t)  ■  ^  Pg  |  Xg,x#,  t  -  (As  cos  a)/vj  (2.1) 

where  v  is  the  medium  velocity  and  a  is  the  angle  of  the  plane  wave  with  the 
horizontal.  The  symbol  As  refers  to  the  spacing  between  sources  and  the 
factor  (As  cos  a)/v  gives  the  appropriate  tiae  delay  between  the  successive 
source  initiation  times.  The  quantity  v/cos  a  is  thus  the  horizontal 
velocity  of  the  plane  wave  and  (As  cos  a)/v  is  the  tine  required  for  this 
wave  front  to  travel  between  two  separate  source  locations. 

The  Fourier-transform  equivalent  of  (2.1)  is 


That  the  imitioB  in  (2.2)  produces  a  plane-wave  sei  snog  ran  nay  be 
seen  by  regarding  the  sanation  as  an  integral  over  the  s  variable  and 
evaluating  that  integral  by  one-dinensional  stationary  phase.  It  is  useful 
here,  however,  to  f  omul  ate  a  nore  comprehensive  approach  in  which  the 
plane-wave  sunation  and  the  nigration  sanation  are  considered  in  unison. 
All  resulting  integrals  can  then  be  evaluated  asynptotically  by  the  nethod 
of  multi-dimensional  stationary  phase. 


Ve  thus  construct  an  inage  or  nigrated  field  in  the  form 

•  CO  00 

Pjjl  (Xq )  -  J  dm  J  d*s  |  dx  Ajj(x#,x  •«»*)  (2.3) 

This  formula  involves  an  integral  over  receivers  followed  by  an  integral 
over  sources.  The  foxier  constructs  the  equivalent  here  of  downward 
continuation  of  the  upward  scattered  wavefronts,  whereas  the  later  provides 
the  plane-wave  construction.  The  integration  over  u  constitutes  the  usual 
imaging  step  of  migration,  where  the  sero-time  data  at  each  depth  level  or 
luge  point  constructs  the  final  migrated  image.  The  vector  Ip  *  (xd,0,zq) 
represents  the  point  of  downward  continuation  where  the  output  image  is 
evaluated. 


From  travel  time  considerations,  the  operator  Ag  is  taken  to  have  the 

fora 


-  4  - 


(2.4) 


Ag  -  A(k.o)  e'ik  [  (ID  '  V  CO,“  +  ZD  ,in  a  +  rD  ] 
with  the  wevenuaber  k  *  w/v. 

Here,  rD  represent*  distance  along  the  path  of  downward  eontinnation. 
The  structure  of  the  phase  in  this  expression  is  based  on  soae  hindsight, 
and  is  constructed  in  order  to  cancel  the  phase  in  the  point  source 
responses . 


The  stationaritjr  conditions  to  be  described  below  will  show  that  the 
exponential  field  of  (2.3)  represents  the  total  phase  of  a  plane-wave  at 
angle  a,  with  the  line  of  zero  phase  passing  through  the  source  point.  The 
xf  cos  a  tens  in  the  exponential  arguaent  gives  the  appropriate  delay  across 
sources  for  the  plane-wave  construction. 


I 

The  result  of  the  stationary  phase  analysis,  the  details  of  which  are 
ontlined  in  the  next  ssction,  is  that  the  aaplitude  factor  A(k, u)  is  chosen 


A(k,  a) 


4«|kjzD  |V  rD  . 

—  1  —  -  sin  a 

rD  rD 


(2.4) 


I  The  geometrical  quantities  involved  here  are  shown  in  Figure  (2).  The 

geometrical  distance  r(  at  first  appears  to  be  unknown  since  it  is  distance 


from  soae  source  location  to  a  specular  reflection  point.  Hence,  it  would 
appear  to  depend  on  a  priori  knowledge  of  the  refleetor.  However,  the 


stationarity  conditiona  illustrated  in  Figure  (3)  show  that  ig  nay  be 
replaced  by 


r  *  (x  -  x  )  coaa+  x_  sin  a  (2.6) 

s  D  a  D 

This  gives  the  correct  geoaetrical  equivalent  for  xg  when  the  iaage  point  Xq 
coincides  with  the  specnlar  point  x*. 

The  above  choice  of  migration  amplitude  gives  an  ontpnt  image  that  has 
a  particularly  simple  fora.  Setting 


r 

g 


=  rD  +  e 


then  as  rD  ->  rg,  s  ->  0.  Evaluation  of  the  spatial  integrals  in  (2.2)  by 
stationary  phase  then  gives 


2jtR  S(w)  e 


i(»/v)e 


+  0(e) 


(2.7) 


Here  S(»)  is  the  Fourier  transform  of  the  source  wavelet,  and  R  is  the 
geometrical  optics  reflection  coefficient.  The  integral  over  <■>  in  (2.2)  can 
be  calculated  exactly,  giving 


PlM(xD>  =  R  S  (e/v) 


(2.8) 


This  lsst  equation  is  the  key  result  of  this  analysis,  and  illustrates 
the  mechanisa  by  which  the  reflector  image  is  created.  Across  any  zone  of 
specular  reflection,  as  shown  in  Figure  (4),  a  family  of  functions  having 
the  fora  of  (2.7)  are  traced  across  the  surface,  with  e  representing  noraal 
distance  to  the  reflector.  The  migrated  output  is  thus  directly 
proportional  to  the  local  geometrical  optics  reflection  coefficient  and  to 
the  area  under  the  Fourier  transform  of  the  wavelet. 

From  the  geometrical  optics  reflection  coefficient,  it  is  possible  to 
infer  some  information  concerning  the  jump  in  velocity  across  ][.  For  plane 
waves  of  incident  angle  0j, 


R 


cos  9^ 


cos 


cos  0T  -  b_ 

1  V 


,  «»»’«!  - 


<2.*> 


Here  &v  is  the  relative  variation  across  Z  of  the  square  of  the  propagation 
velocity: 


6 

v 


(2.10) 


Equation  (2.8)  can  be  solved  for  the  variation  6V  giving 


(2.11) 


4oosaexR  4*in*[(a4$4)/2)R 

(1+R) *  (1  +  R>* 


It  is  then  straight  forward  to  calculate  ▼  in  terns  of  v0  and  6y  by  using 
(2.10).  The  angle  0j,  however,  can  only  be  inferred  from  a  knowledge  of  the 
reflector  tilt,  which  must  be  determined  from  the  imaged  section. 


3.  The  Stationary  Phaaa  Analysis 


In  Appendix  A,  the  construction  of  an  integral  formula  for  the 
reflection  response  at  xg  due  to  a  point  source  at  x$  i*  outlined.  That 
analysis  is  based  on  the  use  of  a  Green's  identity  to  construct  an  integral 
expression  for  the  field  Pg(xg,  xg.<i>)  .  A  Kirchhoff  approximation  is  then 
assumed  for  the  unknown  fields  appearing  in  the  integrand  of  that  result. 
The  integral  over  the  surface  in  the  direction  of  symmetry  of  the  reflector 
is  eliminated  by  one-dimensional  stationary  phase.  leading  to  the 
expression. 


I.*/*  -i  (ff/4)  sgn(k) 

Pc  (x  . x  ,k)  =  -l£J - -j-  e  SU) 

S  *  s  4(2n) 1  * 

.  (3.1) 

N-(*+f>  ik(r  +  r  ) 

R(cr)  .  ■  * - e  8  do 

Jr  r  (r  +  r  ) 

1  s  g  s  g 

7> 


The  integral  here  is  an  integral  in  arc-length  a  across  a  typical  cross- 
sectional  profile  T '  of  the  reflection.  A  combination  of  (3.1)  with  (2.2) 
and  (2.3)  leads  to  a  three-dimensional  integral  in  the  variables  (xg,Xg><r)> 
followed  by  the  integration  over  w: 


PIM(-D> 

OD 

» 

d<u  S '(<*>) 

00 

dz 

8 

00 

dz 

8 

• 

, 

, 

, 

N-(*  +  f  )  ik«(x _.x  . 


(r-Vr. 4 V 


g  s 


where 


i(n/4)  sgn(k) 

S'(w)  =-^J - TT -  S(u>) 

4(2rt)*/ 


and  the  phase  function  f(x^,xs,o) 


is 


€(***»«)  “  r  +  r  -(x.,  -  x  )  cosa-z_  sina  -  r_ 
gs  sgus  D  D 


The  explicit  dependence  of  this  function  on  the  variables  xg, 
given  through  the  distance  functions 


r 

s 


x 


r 

g 


a  (o) 


rD 


where  x  (a)  is  a  point  of  integration  on  the  scattering  surface. 


The  conditions  for  stationarity  in  (xg,x#,o)  of  this  phase 


o)  do 

(3.2) 


(3.3) 

x#  and  o  is 


function  are 


-  10  - 


(3.4) 


i  -4-  (r  +  r  )  =  T-<#  +  !  )  =  0 

a  da  s  g  s  g 


fs  =JT  (rs  +  xsco,a)  = 


-  cos0t  +  cos  a  =  0 


0 

f  =  •* —  (r  -rn)  =  cos  0  -  cos  0„  =  0 

g  3x^  g  D  g  D 

Figures  (2)  and  (3)  show  the  geometrical  quantities  before  and  after 
the  conditions  for  stationarity  are  imposed.  The  first  of  these  is  a 

A 

requirement  for  specular  reflections  (The  symbol  T  refers  to  the  unit 
tangent  vector  along  i’.)  The  second  condition  specifies  that  the  ray  angle 
(with  the  horizontal)  at  the  source  point  must  equal  the  plane-wave  angle  a. 
The  last  equation  is  a  statement  that  the  downward  continuation  angle  0j) 
equals  the  receiver  ray  angle  Og* 


Mathematically,  these  equations  are  specifying  geometrical 
configurations  where  the  dominant  asymptotic  contributions  to  the  three-fold 
spatial  integral  in  (2.2)  occur.  Physically,  the  stationary  directions  are 
the  directions  of  spatial  coherency.  In  the  directions  specified  by  these 
equations,  wavefronts  combine  constructively  in  phase  to  give  a  leading- 
order  contribution.  In  other  directions,  wavefronts  combine  destructively 
to  give  a  lower-order  contribution.  It  is  this  constructive  and  destructive 
interference  that  provides  the  creation  of  (1)  a  net  incident  wavefront 
which  is  planar,  and  (2)  a  backward  extrapolation  of  the  scattered  waves 
which  is  in  exact  reverse  to  the  directions  of  forward  propagation. 


-  11 


The  stationary  phase  formula,  applied  to  a  three-fold  integral  of  the 


fora 


I<X) 


•iX*(5)dx 


is 


I(X)  ~  J 

j 


•  /a 


f{ij) 


i|l|(lj)  +  n/4  (sgn  X)  sig(V^)j 
|det (V j )  | 


(3.5) 


Here  Wj  it  the  aatrix  of  second  derivatives 


»j  *  ■  *  »•*•*  • 

Also,  det(W)  is  the  detexainant  of  this  aatrix  and  sig(Wj)  it  its  signature. 
The  summation  is  over  all  possible  stationary  points  and  the  subscript  on  Vj 
indicates  an  evaluation  at  x  j . 

For  the  problea  at  hand,  the  coapntation  of  the  quantities  det(Vj)  and 
sig(Vj)  can  be  simplified  since,  by  soae  hindsight,  we  anticipate  that  these 
quantities  will  be  of  interest  where 

r  ->  r 

•  I  • 

That  is,  ve  are  interested  in  the  behavior  of  the  spatial  integrals  as  the 
image  point  Xp  of  downward  continuation  approaches  a  specular  point  x '  on 


the  reflector.  The  detail*  of  this  calculation  are  given  in  Appendix  B. 
There,  it  is  shown  that 


|det(W) | 


Zp  sina  cos©j 

~2  *7* 

r  r 
8  * 


(3.6) 


and 


Sig(W)  =  +  1  . 


The  factor  of  cosSj  provided  by  the  presence  of  |det(W)|  in  the 
denominator  cancels  (except  for  a  factor  of  2)  with  ft*(9+9g)  in  the 
numerator  of  (3.2),  since,  at  stationary, 

S>(f  +9  )  *  2cotB  . 

*  8  I 

Also,  the  phase  delay  of  n/4  exactly  cancels  with  the  phase  n/4  in  (2.1). 
The  net  result  of  the  stationary  phase  analysis  is  that 


PIM< 


V  =  } 


A(k,a)R  frg] 
2(k)  sina 


lr  +r 
«  8 


e**  ®  S(«)  dm  . 


(3.7) 


The  explicit  dependence  of  these  variables  is  not  indicated?  however,  it  is 
understood  that  all  quantities  satisfy  the  conditions  for  stationarity  and 
that  (3.7)  strictly  applies  in  the  neighborhood  of  an  image  point.  Also, 
the  summation  over  all  possible  stationary  points  is  understood. 


The  similarity  of  (3.7)  to  a  Fourier  transform  suggests  that  the 
amplitude  factor  A(k,a)  should  be  selected  so  as  to  cancel  the  unwanted 
geometrical  factors  in  that  equation.  Since  onr  interest  is  in  the 


neighborhood  of  an  image  point  on  the  target  reflector,  we  set  r^  =  r^  and 


A(k,a)  =  An  sin  a  jpj 

F? 

(3.8) 

r  +r„ 
s  D 

-  4n  sin  a  sin  0p  ( 

rD 

Also,  Figure  (3)  shows  that  we  may  replace  r(  in  this  formula  by 
r^  -  cos  a-*-  Sin  a 

This  choice  for  A(k,a)  leads  to  the  results  suggested  in  the  previous 
section. 

In  summary,  the  direction  of  this  analysis  has  been  to  assume  a  certain 
canonical  form  for  an  imaging  operator.  This  canonical  form  is  not  specific 
except  for  certain  identified  phase  factors  which  provide  a  reverse 
extrapolation  of  scattered  wavefields.  By  stationary  phase  analysis  the 
amplitude  of  the  image  field  can  be  calculated,  and  this  leads  to  a 
selection  of  the  correct  amplitude  for  the  migration  operator.  This  process 
is  inversion  in  the  sense  that  the  emphasis  is  placed  upon  achieving  a  final 
image  which  depends  only  upon  a  velocity  contrast  or  some  reflection 


coefficient.  However,  it  is  an  indirect  process  in  that  it  does  not  involve 


an  analytical  solution  of  some  scattering  equation. 

The  philosophy  of  the  analysis  involved  here  should  have  useful 
applications  for  other  types  of  exploration  problems.  Specifically,  it  may 
be  useful  for  constructing  migration  operators  which  provide  specified 
cross-sectional  images  from  three-dimensional  data  volumes.  These  operators 
can  be  designed  with  more  consideration  given  to  the  construction  of  true 
amplitude  results.  This  approach  does  not  reqnire  the  migration  operator  to 
be  a  backward  solution  of  the  wave  equation.  It  allows  instead  for  a  more 
general  formulation  of  the  imaging  and  inversion  problem. 
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Conat ruction  of  the  Kirehhoff  Iatogral  Ufnuatitioa 


We  describe  here  the  detsils  of  the  analysis  leading  to  equation  (3.1) > 
which  gives  the  scattered  or  npgoing  response  as  a  surface  integral  over  the 
reflector. 

i 

We  define  first  a  total  aconstic  field  which  is  the  sna  of  the  incident 
and  scattered  fields.  It  is  given  as  a  solution  to  the  equation. 


s  s 


(V  +k  )P  «  -  8(x-Xf)  S(t») ,  k  =  »/v  . 


(A.l) 


The  field  P(x,xg,k)  and  its  noraal  derivative  are  required  to  be 
continnons  across  the  reflector,  and  ?<*,£,, k)  is  also  required  to  satisfy  a 
radiation  condition  at  •: 


liai  R  +  ik  p]  -  0  ,  E  -  |x|  . 


R-»- 


The  scattered  field  Pg(x,5g,k)  is  defined  by  the  equation 


P=PI+PS  • 


where 


ikR 


V*'*.'11  S(*>  •  *,  '  IS-S,  I 


(A. 2) 
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The  function  Pj  is  the  incident  wavefield,  and  represents  the  free-space 
response  to  a  point  source  at  x#.  By  a  standard  Green’s  function  argument, 
Pg  can  be  written  as  an  integral  over  the  surface  of  the  reflector: 


V 


wk)  ■  I  ( ['.<*••  V 


3G(x';x  ) 

k)  — — rs 

dn' 


(A. 3) 


ap,(x',x,,k)  ] 

-5P-  J 


-  G(x'JSg) 


The  free-space  Green's  function  used  here  is 
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and  ft  ia  the  upward  normal  from  the  reflector. 

We  now  invoke  a  Kirchhoff  approximation  in  order  to  simplify  the 
integrand  of  (A. 2).  An  "illuminated”  region  L  on  the  surface  I  is  defined 
as  that  portion  of  the  surface  intersected  by  straight  line  rays  from  the 
source.  The  remainder  of  the  surface  is  denoted  by  D.  We  make  the 


approximation 
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■gjr  <£»£,»*)  =  P(x,x#,k)  *  0  on 


Those  are  the  relationships  predicted  by  geometrical  optics.  The 
reflection  coefficient  R  will  be  taken  as  the  plane-wawe  reflection 
coefficient,  which  for  waves  of  incident  angle  0j  is  given  by  (2.8).  A  sore 
detailed  description  of  the  Kirchhoff  approximation  is  given  by  Bleistein 


The  combination  of  (A. 2)  -  (A. 5)  leads  to 
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Ve  eliminate  the  component  of  integration  in  the  direction  y’  of  symmetry  by 
one-dimensional  stationary  phase  [S],  The  result  is 
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The  final  step  in  the  sinplif ication  of  (A. 6)  is  the  calculation  of  the 
nornal  derivative  of  the  integrand  according  to 


-r  =  V  •* 


Only  the  tern  arising  from  differentiation  of  the  phase  in  the  exponential 
is  kept ,  since  terns  arising  fron  derivatives  of  the  anplitnde  are  of  lower 
order  in  k.  A  farther  precaution  when  sinplifying  (A. 6)  is  to  note  that 
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This  arises  because  a  ninus  sign  has  been  used  in  the  definition  of  the 
forward  Fourier  transforn.  Ve  use  these  results  in  the  stationary  phase 


fornula  (3.S)  and  obtain  (3.1). 


"Calculation  of  tk*  Matrix 


Detfnluat  sad  Sigaature” 


For  application  of  formula  (3.5),  the  determinant  and  signature  of  the 
matrix  of  second  derivatives  most  be  determined.  These  quantities  are 
needed  whenever  the  conditions  for  stationarity  are  met,  hence  equations 
(3.4)  are  applied  in  these  calculations.  This  leads  to  the  following 
quantities: 
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In  the  first  of  these  equations  K  is  the  cnrvatnre  vector  of  X'  at  the 
stationary  point;  K  is  its  magnitude  and  q  =  -  sgn  K*8,  which  is  +1  for  an 
anticlinal  cnrve  and  -1  for  a  synclinal  curve. 

In  addition,  these  quantities  are  needed  as  r^  ->  rg i  i.e.,  near  a 
specular  inage  point  on  the  reflector.  Hence  set  rp  =  rg  and  arrive  at  the 
following  expression  for  the  matrix  of  second  derivatives. 
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The  determinant  of  this  matrix  is  therefore 
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where  cos  0»»  =  Zn/rn 


The  characteristic  polynomial  for  W  is 
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The  following  argument  is  made  for  the  determination  of  the  signs  of 
the  roots  of  this  characteristic  equation.  The  quantity  -ac  gives  the 
product  of  the  three  rootsj  this  quantity  is  always  negative.  Hence,  either 
all  three  roots  are  negative,  or  one  root  is  negative  and  the  remaining  two 
are  positive.  We  then  have  the  following  sub-cases: 


1)  T)  <  0  (reflector  concave  down).  In  this  case  a  +  b  >  0.  But  a  +  b 
gives  the  sum  of  all  three  roots.  This  implies  that  all  three  roots  cannot 
be  negative.  Therefore,  there  must  be  two  positive  roots  and  one  negative 
root,  which  implies  a  signature  of  +1. 
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Figure  (1):  Illustrating  the  reflection  geoaetry. 

Figure  (2):  Illustrating  geoaetrical  quantities  before  stationary 

phase . 

Figure  (3):  Illustrating  geoaetrical  quantities  after  stationary 

phase.  The  angle  at  the  source  is  a,  and  aake  equal  angles  with 
the  surface  noraall  ft,  and  9g  -  G  p  • 

Figure  (4)  Showing  the  iaage  reconstruction  of  a  none  of  specular 
reflection  on  the  target  surface. 


